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Possible formulations of the problem of cosmic rays acceleration in the interstellar 
galactic medium are considered with the use of fractional differential equations. The 
' applied technique has been physically justified. A Fermi result has been generalized 

to the case of the acceleration of particles in shock waves in the supernovae remnants 
fractally distributed in the Galaxy. 
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1 Introduction 
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c/3 ■ The physical meaning of the fractional differential equation 

! 

oDf + D(-A)"/2l^(x,t)=5(x)5;3(t), 6p{t)^t^^/T{l-P), ae(0,2], /3e(0,l], (1) 

ly-^ ! was discussed in previous work [1] in the continuous time random walk (CTRW) model in 

OO I view of the problems arising when it is applied to the description of the transport of cosmic 

^ I rays in the Galaxy. In this work, the CTRW approach is considered in application to the 

• description of the acceleration (more precisely, reacceleration) of cosmic rays. In contrast 

O ' to the preceding problem, this concerns the behavior of particles in the momentum 

space. Successive interactions (collisions) of a charged particle with more or less localized 
inhomogeneities of the magnetic field from slowly moving magnetic clouds mentioned 
by Fermi in his pioneering works [2] to strong shock waves in the supernovae remnants 
rS ■ mentioned by Berezhko and Krymskii in review [3] can be considered as instantaneous 

■ jumps from one point of the momentum space to another one. The momenta Apj acquired 

by the particle in these collisions are random and, even for their isotropic distribution, 
the point 

p = Po + Api + Ap2 + Ap3 + . . . , 

representing the particle in the momentum space moves away from the point (momentum) 
of the acceleration injection po similar to a Brownian particle; this behavior means the 
further acceleration (reacceleration) of the particle. However, only a certain fraction of 
the particles moving away from the center are accelerated. This fluctuation component of 
the mechanism of the acceleration of cosmic rays is analyzed in this work. 

From the statistical point of view, the main consequence of the Fermi conclusion is that 
the exponential increase in the energy of the accelerated particle with time, E = Eqc""^, 
and the exponential distribution of the age of the detected particles, dP = exp(— t/r)(it/r. 
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are sufficient for tlie formation of the power-law energy spectrum N(E). Tliat is all. There 
are no other sources of fluctuations taken into account in the Fermi model. This result 
is evident: What fluctuations can else exist if an increase in the energy by a factor of e 
requires 10^ coUisions according to the Fermi estimation? A more significant effect could 
be produced by more rare acts with more large acceleration in each of them. Among these 
processes are the aforementioned interactions with strong shock waves, when even a single 
interaction can increase the energy of the particle by a factor of 7 -^ 13 (see [4, p. 449]). 
To this end, it is appropriate to pass from the degenerate spectral function d{E — Eqc"'^) 
characterizing the determinate process of the Fermi acceleration to the continuous function 
n{E, t) related to the momentum distribution /(p, t) as 

n{E,t)^ J S{E-E{p))f{p,t)dp. 

Under the Fermi assumption that the parameter r is independent of the energy, the desired 
spectrum can be represented in the form 



N{E) = N{E;t)^t- 

where 



J n{E,t)e-*/^dt 

-0 



r-'n{E,T-'), (3) 



n{E,\) = J e-^^n{E,t)dt 



is the Laplace transform of the spectral function in the time variable. The effect of 
the kinetic fluctuations is now taken into account at the stage of the construction of 
equations for the distributions /(p, t) or n{E,t) by the inclusion of additional terms 
containing differential and integral operators. The theory of the acceleration of cosmic 
rays is under development and this work is focused on the role of fractional derivatives in 
this development. 



2 Fractional kinetic equations 

As was mentioned in [1], Eq. (1) is the diffusion limit of the coordinate (x — CTRW) 
model, where the trajectory of a particle is a sequence of instantaneous jumps at random 
times by random distances; the particle is at rest between these times. The CTRW model 
with a finite velocity of free motion provides qualitatively new results (see, e.g. [5]). 
When analyzing the acceleration of cosmic rays in the first approximation, the stochastic 
dynamics of the particle in the velocity (momentum or energy) space is considered, and 
the boundedness of the coordinate space is taken into account by correcting the mean age 
T or introducing additional coefficients characterizing the leakage of the particles from the 
acceleration region or from the Galaxy as a whole. 

In the momentum space, the stationary position of the point that does not coincide 
with the coordinate origin means the motion of the particle with a constant momentum 
(velocity, energy). The exponent /3 characterizes now the tail part of the distribution of 
the random interval duration AT between the successive collisions of the moving particle 
Q{t) = P(AT > t), i — >■ oo. For the ultrarelativistic particle {v ~ c), this interval is 
proportional to its mean free path between the successive collisions, so that the new (3 
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values should correspond to the a values in the preceding x — CTRW model. In view of 
this circumstance, the consideration starts with the CTRW system of Eqs. (2)-(4) from 
[1], where the coordinates r are changed by the momenta p, the momentum change rate 
V is accepted to be infinite, and it is assumed that the particle is injected by the source 
at the time i = with the distribution function /o(p): 

t 

/(p, t)^ j Q{t - t')Fo^r(p, t')dt', (4a) 



Fo^i(p,t) = Jw{p^ p')F,^o{p',t)dp' + Mp)S{t), (46) 

t 

Fi^o(p, t) = J q{t - t')Fo^i{p, t')dt'. (4c) 



The further transformation of the system is associated with the specification of the 
distributions q{t) and w{p ■<— p') distributions. The time is usually (one can say, always) 
taken in the exponential form Q{t) = e~^*,g(t) = —dQ/dt = ^e~^^ and therefore the 
master equation begins with the first-order time derivative df{p,t)/dt = . . .. This means 
that the process is assumed to be Markovian. However, the real distribution of the time 
intervals between collisions is unknown. For example, one can assume that it is of power 
law, Q{t) oc t"", rather than exponential. This is in agreement with the self-similar pattern 
of turbulent motions and with its power-type laws. The hypothesis of the fractal character 
of the interstellar medium [6] also provides power-law distributions. This concerns the 
behavior of magnetic field lines in interstellar space. They are usually represented as 
relatively smooth lines, which somewhere rest at the magnetic traps, intersect each other, 
sharply change their directions, and performing a "diffusion dance"in time. The leading 
centers of particles moving along spiral trajectories move along these lines. If these smooth 
sections become invisible and chaotic patterns of the structure become prevalent with the 
expansion of the field of view and the corresponding decrease in the scale, then this is 
an asymptotically homogeneous medium where the mean free paths can be simulated by 
a usual exponential function. If the expanding field of view includes increasing straight 
segments replacing the segments becoming small due to a decrease in the scale so that 
the structure remains (qualitatively) unchanged, this is a fractal structure. Under these 
conditions, the exponential distribution of mean free paths characteristic of a strongly 
mixing medium cannot be expected, but can not be completely rejected as well. The best 
compromise would be a family of distributions including both exponential and power-law 
distributions. Fortunately, such a family exists: it is a set of the functions 

Qa{t) = E^{-lj,r), aG(0,l], 

oo 

where Ea(z) — /V{an -\- 1) are the Mittag-LefSer functions. The function Qa{t) 

n=0 

with q; = 1 is a usual exponential and with a < 1 is a fractional exponential having a 
power-law asymptotic behavior t — )■ cxo. The corresponding density qa{t) satisfies a 
fractional differential equation [7]; as a result, the entire system of Eqs. (4a)-(4c) has the 
fractional time-differential form 

oDr/(p, t) = //A/(p, t) + h{p)5^{t). (5) 
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Here 

A/(p, t)^l w{p ^ p')/(p', t)dp' - /(p, t) (6) 

is the acceleration integral with transitions density ^(p ■<— p'). The time series of 

accelerating collisions forms a fractional Poisson process of order a [7], with a — )■ 1 
becoming an ordinary Poisson process which underlies classical kinetic equation (see 
(21.1)in [8]). Investigations performed in [9] indicate a qualitatively new property of this 
process: the average number of collisions increases more slowly (oc t") than in the usual 
case (oc t) and relative fluctuations of the number of collisions in the limit t — > oo do not 
disappear, but tend to a limiting distribution depending on a (some sort of KNO scaling). 

The transition from the first time derivative to the fractional derivative of the order 
a < 1 does not require solving a fractional differential equation; it is more convenient to 
use the relation between the solutions of Eq. (5) with a < 1 and a — 1 [7], 

oo 

/„(p,t) = {t/a) J /i(p,r)<7+(tr-i/";a)T-V"-MT. (7) 



Here g+{t;a) is the one-sided stable (in the Levy sense) probability density determined 
by the Laplace transform 

/■oo 

/ exp(—Xt)g+(t; a)dt — exp(— A"), 
Jo 

and t, T , A are dimensionless variables. The Laplace transform of Eq. (7) in time with 
the use of the above formulas for the spectra provides the formula 

N^{E-t)^N,{E-t"). (8) 

It presents the effect of the fractal dimension of the fractional Poisson process of collisions 
a e (0, 1] on the energy spectrum of cosmic rays: the spectrum Na{E;T), formed by an 
ensemble of particles with the mean life time r, which are accelerated according to the 
fractional Poisson law of the order a < 1, coincides with the spectrum of particles that 
are accelerated by the usual Poisson process (a = 1), but have the mean age r (recall 
that the time is here dimensionless and the injection spectra fo{E) are the same in both 
problems). On the example of the Fermi spectrum, it is easy to see that the efficiency of 
acceleration decreases (the spectrum becomes steeper) with a decrease in the order of the 
process. The fractal character of the spatial distribution of accelerating regions naturally 
reduces the efficiency of acceleration. 

3 Fractional Fokker-Planck equations 

Similar to the classical case, the transition from the kinetic equation to the Fokker-Planck 
equation is associated with the transformation of the collision integral to a differential 
form by expanding the integrand into a series in the momentum increment to the second 
order terms. There are two variants of such an expansion, which provide slightly different 
equations (see, e.g., [8]). The first variant implies the smallness of the absolute value of the 
change in the momentum | Ap| = |p — p'|, so that the momentum of the incident particle 
only slightly changes in the magnitude and direction in a single collision event (e.g., as in 
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the case of the coUision of a heavy particle with a hght one). The second variant imphes 
the smallness of the change in the absolute value of the momentum A|p| = |p| — |p'|, 
whereas the change in its direction is not small and can have a wide distribution up to 
the isotropic one (in the case of the collision of a hght particle with a heavy one) . Under 
the assumption of the isotropic scattering, the fractional differential generalization of the 
Fokker- Planck equation is obtained in the form 

oD^ f{p,t) = Ap{K{p)f{p,t)) + Mp)5a{t), (9) 

where 



K{p) = {tx/2) J {Apfw{Ap; p)d Ap 



is the diffusivity in the momentum space. An energy analog of Eq. (9) {a — 1) is well 
known in the physics of cosmic rays in the form (see Eq. (14.2) in [10]) 

r.a,(p ,. d[a,{E)n{E,t)] d'[a2{E)n{E,t)] 
o^t t) \ \- no[E)da[t). (lUj 

At the same time, Eq. (9) with diffusion term 

Ap(i^(p)/(p, t)) = (Api^(p))/(p, t) + 2(VpX(p)) Vp/(p, t) + K{p)A^f{p, t) 

significantly differs from another diffusion type equation (cm. (9.57), [10]) 

oDt f{p,t) = Vp(i^(p)Vp/(p,t)) + fo{p)Sa{t). (11) 

The difference is due to the fact that Eq. (9) is derived in frame of the coUision model, when 
the point presenting the particle instantaneously moves to another, possibly far, geometric 
point, violating the continuity of the trajectory in the momentum space, whereas the 
dynamic derivation of Eq. (11) implies that the trajectory in the momentum space is 
continuous and even differentiable. The classical versions (with a = 1 ) of Eqs. (9)-(ll) 
underlie the standard set of mathematical tools describing the fluctuation mechanisms 
of the acceleration of cosmic rays; their solutions are well known [4, 10, 11]. A common 
property of these models is the Gaussian character of the momentum distributions, which 
is due to the assumption that the second moment of the momentum acquired in an 
acceleration event, which enters into the diffusion coefficient, is finite and the acceleration 
rate is low. A change of the first time derivative in this equation to its fractional analog 
of the order a e (0, 1) does not lead to an increase in the efficiency of acceleration; on 
the contrary, the acceleration rate in the subdiffusion regime decreases further. It is also 
noteworthy that the statement of some authors that diffusion at a > 1 is accelerated is 
erroneous. This statement is based on linguistic intuition ("if a < 1 means subdiffusion 
and a — 1 corresponds normal diffusion, then a > 1 should mean super diffusion") and 
is erroneous, because the solution f{p,t) at a > 1 loses its probability meaning: it is 
not positively definite in this case. Thus, the transition to the fractional time derivative 
(with the unchanged remaining, diffusion part of the equation) as a method for enhancing 
the high energy part of the spectrum is physically unpromising. At the same time, the 
parameter a presents the existence of possible correlations in the spatial distribution of 
the acceleration regions and can be useful in this respect (recall that a — 1 corresponds to 
the uniform Poisson distribution of such regions, which do not correlate with each other). 
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This parameter can be kept and the efficiency of the multiple acceleration mechanism 
can be increased only by modifying another operator of the equation, more precisely, 
by returning from the differential form of the acceleration operator, which describes a 
continuous slow collection of the energy, to the integral form describing acceleration as 
a sequence of events with large instantaneous (in terms of the considered "galactic "time 
scales) changes in the momenta. For this reason, it seems appropriate to change the 
momentum Laplacian to its fractional analog, because the fractional Laplacian contains 
a momentum integral operator with inverse power law kernel, which can ensure a high 
acceleration rate. Let us consider this variant in more detail. 

The losses of the energy of a fast charged particle in a medium, which is described 
by equations similar to Eq. (10) (naturally, with a — 1), cannot be larger than its 
initial energy and all of the moments of the energy loss are finite. In the problem 
of acceleration, there is no such definite limit of the energy increase; this fact is an 
additional reason to study the region with infinite dispersion, which attracts an increasing 
interest of researchers of anomalous diffusion processes. In this case, informative (in the 
asymptotic sense) results are obtained only when infinite dispersion is due to the power- 
law distributions: 

J w(Ap; p')rfAp oc p"^, p ^- oo. (12) 

Ap|>p 

If 7 > 2, the second moment is finite and this corresponds to the classical diffusion 
region. If 7 < 2, the second moment of the increment is infinite and this corresponds 
to the model of additive Levy fiights. In this case, the equations for the momentum and 
energy distributions f(p, t) and n(E, t), respectively, following from the asymptotic CTRW 
analysis, 

oDffip,t) = -Ki-A^y/'f{p,t) + Mp)6^it), (13) 

and 

r)a^(p f^ _ / d^ME,t)]/dE- + no{E)5^it), < u < 1; 

° * ' ^'~\ d[ain{E, t)]/dE + a"[a,n(E, t)]/dE'' + no{E)Sa{t), l<v<2. ^'-^> 

seem to be reasonable. Here 

; • ; ^ ^! 



V — 



2, 7>2, 



jdE^ is the symbol closer to the ordinary notation of the fractional differential operator 
and a\ and are the constant coefficients. The constancy of these coefficients is very 
significant for the derivation of the equations. Fractional differential equations are usually 
derived with the use of integral transformations whose applicability requires the constancy 
of the coefficients. It would be incorrect to derive, e.g., Eq. (13) in such a manner and, 
then, to place the variable diffusion coefficient K(p) in front of the fractional Laplacian 
(this is obvious even on the example of Eqs. (9)- (11) with the integer Laplacian). 

Equation (14) is clearly similar to the i^th term approximation of the expansion of 
a function that is zero at the reference point and it would seem that omitting the next 
(divergent!) term introduces an infinite error, but this is not the case. The reason is that 
the Taylor formula rather than the infinite Taylor series is used and the series expansion 
can be continued only until derivatives exist; after that, the remainder term, which is 
always finite, should be written. If the next derivative existed, its inclusion in the continued 
expansion would mean an approximation of the omitted term. If this derivative does not 
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exist, it cannot be used to approximate the omitted term and one should return to the 
initial point, where the derivative of the preceding order is used for approximation. The 
diffusion packet, which is described by Eq. (13) and propagates from the origin of the 
momentum coordinates, has the form of the three-dimensional isotropic fractional stable 
distribution G (0,2], aE (0, 1], smearing proportionally to t"''": 

fip,t) = {Krr'/^4'^'''\{Kn-'/^p). 

Tails (or, more romantically, wings) of this distribution have a power-law form with the 
exponent u. Physically, this means a peculiar leading effect: one of independent terms in 
the sum Api + Ap2 + . . . + Ap„ is always outstanding and compared in magnitude with 
the reminder. This leading effect disappears at i> = 2, when the distribution becomes 
Gaussian (sub-Gaussian). As a result, the spectrum at u < 2 has the form 

Ni{E)dE (X E-'^-^dE, 

which is similar to the Fermi formula. The main difference is that the exponent u in this 
case is independent of the age of the detected particles and is completely determined by 
the acceleration mechanism in an individual local event (collision). For this reason, the 
fractal character of the spatial distribution of accelerating regions also does not affect the 
slope of the resulting spectrum. 

4 Integro-fractionally-differential model 

A drawback of the model with the fractional momentum Laplacian is that the increments 
of the momentum in the acceleration event are independent of the momentum of the 
particle involved in the interaction, whereas in the Fermi model and in its later variants, 
the increments of the energy (and, therefore, the momentum) are on average proportional 
to the energy (momentum) of the particle before the interaction. In this case, the energy 
of the accelerated particle is expressed in terms of the product of independent random 
variables, rather than their sum. This model is called multiplicative walk in order to 
distinguish it from the additive walk model considered above. The increment of the 
momentum in the multiplicative model is proportional (in the statistical case) to the 
absolute value of the momentum p' of the particle coming into interaction, 

Ap = p'q, J w{Ap; p')dAp oc {p/p')~^, p ^ oo. (15) 

|Ap|>p 

Under the assumption that the distribution of the proportionality vector q is independent 
of p' and isotropic, 

iy(q; pOdq = {l/2)V{q)dqdi, i = cos(q, p) 
, kinetic equation (5) can be modified to the form 

{1 , oo si 
I il + + / (Vl + 2e? + 5') dq-f{p,t)\+U{p)5o.{t), 

(16) 
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representing a new model of a distributed reacceleration, more precisely, a new 
modification of the model proposed in [12, 13]. In order to make this model be closer 
to real processes of reacceleration, e.g., in the case of the intersection of shock fronts in 
the remnants of supernovae, we assume that [13] 

The resulting model can be called multiplicative Levy flights. 

Let us consider the equation for the spectral function in two extreme cases. In the first 
case, 7 > 2, the second moment of the momentum increment proportional to exists, 
and this case corresponds to the classical diffusion with variable coefficients: 

d[a,En{E,t)] d'[a2EME,t)] 
oU^n(ii/,tj — \ \-no{E)da[t). 

In the second case, we suppose 7 << 2, so that only the term may be retained in the 
radicand in Eq. (16): 

oD^iE, t) = /X I Jjq-^-'n{E/q, t)dq/q - n{E, i) | + no{E)5^{t). (17) 

I do not claim that this approximation is very good, but namely this acceleration operator 

was used in [13] for some calculations. Solving Eq. (17) with the method of the Mellin- 
Laplace transforms and using Eqs. (3) and (8), the following expression is obtained for 
the case of the mono-energetic source {n{E) = 6{E — Eq)): 

N^{E;t)= , —) — . (18) 

' ' (1 + /XT")2 \EoJ Eo ^ ' 

Although Eq. (18) was derived under the assumptions that strongly simplify the real 
situation and has a qualitative sense, it compactly presents the effect of all three sources of 
fluctuation acceleration: fluctuations of the age of the particle (parameter r), fluctuations 
of the number of the acceleration events (a h ^), and fluctuations of the energy acquired 
in a single event (7). Representing the scaling parameter 11 in the form = r^, where ta 
is the characteristic time interval between the events of the acceleration of particles in the 
remnants of various supernovae (recall that r is the mean lifetime with respect to nuclear 
collisions), the absolute value of the exponent of the integral spectrum can be written in 
the more clear form 7' = 1 + 7/(1 + (t/ta)"]. At q; = 1, and /ir ^ 1, we arrive at the 
Fermi formula (2) with a = \xj^. 



5 Conclusions 

Taking into account that all three above-considered mechanisms contribute to the 
acceleration (more precisely, reacceleration) of cosmic rays, this process can be analyzed 
using the integro-fractional-differential equation 

oD,"/(p,t) = /xoVp(iro(p)Vp/(p,t)) +/xiAp(iri(p)/(p,t)) -/X2ii'2(-Ap)^/V(p,t)+ 
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I / Y Jv{q)f (p/^Jl + 2iq + q\ / [^^l + 2^5 + g^)' rfg - + h{p)5^{t). 

(19) 

The terms with the coefficients /io, yU-i, /i2 and /i3 successively present the contributions 
from the processes of continuous acceleration (such as fluctuations of the momentum 
of the particle in turbulent plasma [11, 14]), collision accelerations with a finite second 
moment (the model of walking magnetized clouds with a limited distribution of their 
velocities), additive Levy-type collision accelerations given by Eq. (12) (the same with 
the power-law distribution of the velocities of the clouds), and multiplicative Levy-type 
collision accelerations given by Eq. (15). Recall again that the last operator describing, in 
particular, acceleration at shock waves in the remnants of supernovae does not reduced 
to the fractional differential form, but holds in the integral form. Equation (19) should be 
supplemented by terms presenting energy losses, exist of particles from the Galaxy, and 
nuclear interactions (fragmentation). 

However, only dominant terms can be retained in Eq. (19) applied to particular 
problems, as is usually done. 
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